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The development of fractures at regular length scales is a widespread feature of Venusian 
tectonics. Models of lithospheric deformation under extension based on non-Newtonian viscous flow and 
brittle-plastic flow [1 ,2] develop localized failure at preferred wavelengths that depend on lithospheric 
thickness and stratification. The characteristic wavelengths seen in rift zones and tessera can therefore 
provide constraints on crustal and thermal structure [3]. Analytic solutions were obtained by [4] for growth 
rates of infinitesimal perturbations imposed on a one-dimensional, layered rheology. Brittle layers were 
approximated by perfectly-plastic, uniform strength, overlying ductile layers exhibiting thermally-activated 
power-law creep. This study investigates the formation of faults under finite amounts of extension, 
employing a finite-element approach. Our model incorporates non-linear viscous rheology and a Coulomb 
failure envelope in the manner of [5]. An initial perturbation in crustal thickness gives rise to necking 
instabilities. A small amount of velocity weakening [cf. 6] serves to localize deformation into planar regions 
of high strain rate. Such planes are analogous to normal faults seen in terrestrial rift zones. These "faults" 
evolve to low angle under finite extension [cf. 7], Fault spacing, orientation and location, and the depth to 
the brittle-ductile transition, depend in a complex way on lateral variations in crustal thickness. In general, 
we find that multiple wavelengths of deformation can arise from the interaction of crustal and mantle 
lithosphere. 

To calculate patterns of deformation, we use a finite element (FE) program for incompressible 
Newtonian viscous fluid flow via the penalty method [8], modified to treat non-Newtonian flow problems 
[9,10], Gravitational body loads are applied to individual elements. An initial viscosity is calculated for a 
state of uniform extension using the power-law relation 

Heft = Be max ^ ’ n ) /n exp(Q/nRT), (1 ) 

where is the maximum principal shear strain rate, n is the power-law exponent, R is the gas constant, T 
is temperature, Q is the material activation energy, and B is a material strength constant. An iterative 
procedure re-calculates the effective Newtonian viscosity |J.(i,j) for each quadrilateral FE. The maximum 
viscosity is limited to that allowed by the Coulomb failure criterion. Before proceeding with any time-steps, 
we iteratively calculate the deviatoric strain rate at zero strain resulting from the brittle-plastic and velocity- 
weakening effects. Strain rates are locally enhanced in elements with ductile strength less than that of 
mantle at the base of thickened crust. Higher strain rates imply reduced effective viscosity of these 
elements. Strain is preferentially accomodated where viscosity is lower, leading to further reduction in 
strength. This feedback eventually results in highly localized rates of deformation, i.e., the material breaks 
and undergoes faulting. Localization of deformation is the result of non-linear dynamics and is somewhat 
unstable with respect to initial conditions. Generally more than 100 iterations are required to reach steady 
state, even with over-relaxation on strain rate. Finite deformation is then calculated by a second-order 
time-step integration, with dynamic topography developing at the free surface. 

For the crust we use a recent experimentally obtained flow law for the rheology of diabase [1 1 ] 
under extremely dry conditions representative of the surface of Venus, with n=4.3 and Q=482 kJ/mol in 
equation (1). The stronger olivine rheology preferred by [5] is assumed for the mantle, with n=3 and 
Q=540 kJ/mol. The brittle-ductile transition thereby occurs at about 8-10 km in the crust and 15-20 km in 
the mantle. We employ a lithospheric thermal model as discussed in [3], with a surface temperature 
gradient of 25 °K/km. The initial strain rate is Eq = 10‘ 14 s' 1 . A brittle strength evelope provides an upper 
bound to viscosity at a given depth and strain rate, such that deviatoric stress T max ^ o 0 + Cpgz. The 
Coulomb strength coefficient C is velocity-weakening, depending on strain rate as C=C 0 - 0.05 logio(£/£o). 
with 0.25<C<0.4, C 0 = 0.36 and a 0 =1 MPa. This value of C 0 corresponds to an internal friction coefficient 
of about 0.71 . A perturbation of crustal thickness between ±1 and ±10%, at a 50 to 100 km wavelength, is 
introduced. 

For a range of crustal thicknesses of 8 to 18 km, preliminary results show that faulting develops at 
multiple scales when brittle-ductile transitions occur in both the crust and upper mantle, with the width of 
the shortest wavelength features being about twice the thickness of the brittle crust. Longer wavelength 
features occur at twice the depth to the base of the upper mantle lithosphere. The model boundary 
undergoes a prescribed rate of uniform extension. Integration of strain rate over time to overall strains of 
10% produces conjugate patterns of faulting and rifting. Larger amounts of strain generally exceed the 
range of applicability of quadrilateral FE's. 
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The pattern of normal faults and sloping terraces (Figure 1) indicates a fairly realistic simulation of the 
behavior of layered brittle-ductile rheology. Localized faulting develops when the initial perturbation 
wavelength is less than about 8 times the total thickness of the lithosphere (base of brittle mantle). 
Longer wavelength perturbations produce stable extension of the crust, and brittle-plastic necking of the 
mantle lithosphere. Perturbations at wavelengths of less than twice lithospheric thickness produce a 
single pair of conjugate normal faults meeting at the base of the crustal lithosphere. The growth of fault- 
like zones indicates a mechanism for the development of linear features at multiple scales due to 
extension in a Venus lithosphere with a variable-thickness crust. We plan to use this model to investigate 
the relationship between deformation patterns, wavelength of perturbation, and strength envelopes. 
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Figure 1 . Finite element model with assumed flow laws and thermal structure. Original crustal structure is 
50 km wide, 15 km deep, with ±5% variation in crustal thickness, greatest at the center. Density of crust is 
2900 kg/m and mantle is 3300 kg/m 3 . The center of the grid is a symmetry axis, with periodic boundary 
conditions on the sides. The bottom boundary at 50 km depth has no vertical displacement, and is shear- 
stress free. The grid shown has undergone a total extension of 6%. Regions of high strain rate are shown 
on the left by the darkly shaded regions. Regions on the right are shaded where the Coulomb brittle 
failure envelope is locally exceeded. Dynamically produced topography at the surface reflects the 
perturbation wavelength as well as finer scales of faulting generated by the interaction of crustal and 
mantle lithosphere. 


